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ABSTRACT 

Measuring the three-point correlators of the Cosmic Microwave Background (CMB) anisotropies could 



■ help to get a handle on the level of non-Gaussianity present in the observational datasets and therefore 

would strongly constrain models of the early Universe. However, typically, the expected non-Gaussian 
^ ' signal is very small. Therefore, one has to face the problem of extracting it from the noise, in particular 

' from the 'cosmic variance' noise. For this purpose, one has to construct the best unbiased estimators 

for the three-point correlators that are needed for concrete detections of non-Gaussian features. In this 
article, we study this problem for both the CMB third moment and the CMB angular bispectrum. We 

O ! 

emphasize that the knowledge of the best estimator for the former does not permit one to infer the best 
estimator for the latter and vice versa. We present the corresponding best unbiased estimators in both 



cases and compute their corresponding cosmic variances. 



Subject headings: cosmic microwave background — methods: analytical — cosmology: theory — large 
' scale structure of universe — early Universe 

1. INTRODUCTION 

fflhe Cosmic Microwave Background (CMB) has been recognized as one of the best tools for studying the early Universe 
^ • (e.g. Scott 1999 ). In particular, the statistical properties of the CMB anisotropies are a powerful means to discriminate 
amongst the possible scenarios. This is because, in general, different models predict different statistical properties. For 
example, the simplest models of inflation predict that the temperature anisotropies should obey a Gaussian statistics and 
therefore any non- vanishing measurement of a three-point correlator (in a sense to be precised below) would automatically 
ruled out such models, a very interesting result indeed. 

From a practical point of view, measuring any non-Gaussianity in the data is a very difficult task since the signal is 
typically very small. Of course, this signal should be compared to the noise and what really matters is the signal to noise 
ratio. The noise can have many different origins including instrumental errors, foregrounds contamination or incomplete 
sky coverage. Another source of error is the so-called 'cosmic variance'. Roughly speaking, it comes from the fact that we 
only have access to one realization of the temperature anisotropies whereas theoretical predictions are expressed through 
ensemble averages. In a Gaussian model, for example, the mean value of any three-point correlator has to vanish but this 
does not guarantee that a concrete detection of a non-zero signal on the sky would be in contradiction with the model 



(Scaramella & Vittorio 1991, Srednicki 1993). The important point is that the cosmic variance can dominate the other 
sources of error, as this is in fact the case for the two-point correlators on large angular scales. Therefore, if one wants to 

1 A.G. would like to dedicate this article to the memory of his Ph.D. supervisor, Dennis William Sciama. 
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unveil non-Gaussianity, it is necessary to address the cosmic variance problem for the three-point correlators. The usual 
way to deal with this problem is to construct estimators by performing spatial averages on the celestial sphere and to 
find the one which has the smallest possible variance. The aim of this paper is then to find the best unbiased estimators 
both for the third moment (a™ 1 a^ 2 a^ 3 ) and for the angular bispectrum Cg 1 e 2 e 3 and to display the corresponding cosmic 
variances. 

Recently, there has been a lot of activity in the subject triggered by the finding that non-Gaussianities are present in 



the 4-yr COBE-DMR data ( Ferreira et al. 1998 , Pando ct al. 1998 ). Further analyses have confirmed this result (e.g. 



Bromley fc Tegmark 1999). However, soon after, it was demonstrated by Banday et al. (1999) that the non-Gaussian 



signal is driven by the 53 GHz data. This systematic artifact in the CMB maps rejects a possible cosmological origin. 



More generally, it is clear that the presence of foregrounds (Bouchet & Gispert 1999, Tegmark et al. 1999) renders difficult 
the detection of a genuine non-Gaussian signal. Nevertheless, one should expect non-Gaussian features to be present in 
the CMB anisotropy datasets. These could be produced in the early Universe during inflation either because the initial 



conditions are non-Gaussian themselves [i.e. the quantum initial state is not the vacuum (Martin ct al. 1999, Contaldi et 



al. 19^9J)] or owing to the existence of couplings between different perturbation modes at the non- linear level (Gangui et 
al. 19^4. Gangui 1994 , Lindc fc Mukhanov 1997 ). In the context of slow-roll inflation, the CMB bispectrum has recently 
been studied in ( jGangui' fc Martin 200C| , |Wang fc Kamionkowski 2000| ). Even if non-Gaussianities are not primordial 
in origin, they will nevertheless arise during later stages of evolution. In this context, the Rees-Sciama effect will build 
up a small but non-vanishing signal (Luo & Schramm 1993, Mollerach et al. 1995, Munshi et al. 1995). Also, cosmic 
topological defects of the vacuum, like strings and textures, are amongst the best motivated sources for non-Gaussian 



features (Bouchet ct al. 1988, Ferreira fc Magueijo 1997 , Avelino et al. 1998, Gangui fc Perivolaropoulos 1995, Gangui 



Molle f-ach 1996[ ). Regarding secondary sources, Goldberg & Spergel (1999) and Spergel & Goldberg (1999) have recently 
calculated the angular bispectrum due to second order gravitational effects like the correlation of lensing of CMB photons 
and secondary anisotropics coming from the Integrated Sachs Wolfe effect and thermal Sunyaev-Zel'dovich effect. In 
the same line, Cooray and Hu (1999) have taken into account further additional contributions to the bispectrum in the 
presence of reionization. Other approaches to the study of non-Gaussian features include preferred-direction statistics 



for sky maps ( Bunn fc Scott 1999 ), the three-point correlation function ( Falk et al. 1995 , Hinshaw et al. 1994 , Gangu: 
et al. 1994 , Hinshaw ct al. 1995 ), lensing statistics ( Bcrnardcau 1997 , Winitzki 1998 , Zaldarriaga 1999| ), the genus and 
Euler-Poincare statistics (Coles 1989, |Gott ct al. 1990 , Bmoot 1999] ), peak statistics (Bond & Efstathiou 1987, Kogut ct 



al. 19j9q, Kogut ct al. 1996), correlation function of peaks (Heavens & Shcth 1999), Minkowski functionals (Winitzki & 
Koso^ fsky 1998 ) and wavelet analyses ( Popa 1998 , Hobson et al. 1998 ). 

This article is organized as follows. In the next section, the general strategy for finding best estimators is exposed. As 
a warm up, in the third section, we implement this strategy for the two-point correlators. The fourth section is the core 
of the article. There we explicitly derive, for the first time, the best unbiased estimator for the angular bispectrum and 
show its corresponding variance. Except for an overall normalization factor, this estimator turns out to be the one already 
employed by Ferreira et al. (1998) and other authors recently. Our result places their choice on a firm basis. Next, we 



find the expression for the best unbiased estimator for the third moment. An earlier study was performed in (Heavens 
19980 however, our findings go beyond the results obtained in that article and, moreover, are explicit. Moreover, we 
present its corresponding variance. In addition, we also emphasize that the knowledge of the best estimator for the third 
moment does not allow one to infer the best estimator for the angular bispectrum and vice versa. In the last section, we 
briefly present our main conclusions. We finish up with a short Appendix which includes formulae related to the inverse 
two-point correlation function. 



2. GENERAL STRATEGY FOR FINDING THE BEST ESTIMATOR 



In this section, we expose the cosmic variance problem from the viewpoint of the theory of cosmological perturbations of 
quantum-mechanical origin and describe the method of the best unbiased estimators. This theory rests on the principles 



of general relativity and quantum field theory. At the beginning of the inflationary phase ( Guth 1981 , Lindc 1983a 



Albrecht & Stcinhardt 1982, Lindc 1983b) the Friedmann-Lemaitre-Robertson- Walker background spacetime already 



behaves classically whereas the excitations of the metric around this background are still quantum mechanical in nature. 
Technically, this means that the perturbed metric must be considered as a quantum operator. This operator either 
represents density perturbations or gravitational waves. In each case, the quantization can be carried out in a consistent 



way (Mukhanov fc Chibisov 1981|, Hawking 1982, Starobinsky 1982, Bardeen et al. 1983, Mukhanov et al. 1992, Grishchuk 



1993, Martin & Schwarz 1998, 1999). Then, the (zero-point) quantum fluctuations, which are the seeds of the cosmological 



perturbations, are amplified during inflation owing to the particle-creation phenomenon or squeezing effect ( Grishchuk 



& Sid arov 1990 ). Next, these primordial fluctuations give rise to the large scale structures and to the CMB anisotropics 
observed today in our Universe. 

One should also discuss the choice of the quantum state in which the metric operator is placed. Obviously, it is not 
possible to prepare the initial state of the Universe and therefore the choice of the quantum state of the perturbations is 
a priori free unless some theory of the initial conditions is provided [for example, quantum cosmology ( Halliwcll 1989| )]. 



Usually, it is assumed that the initial state is the vacuum although different hypothesis are possible ( Brandenberger & 



Hill 1 )86, Martin et al. 1999, Contaldi et al. 1999). If the initial state is the vacuum, then the corresponding statistical 



properties are Gaussian. This is because the ground-state wave function of an harmonic oscillator is a Gaussian. It is 
possible to avoid this general conclusion either by considering non-linear cosmological perturbations or by assuming that 



the initial state is a non- vacuum state. We have recently investigated the first possibility in (Gangui & Martin 200C). The 
second possibility has been studied by Martin et al. (1999). In the latter case, non-Gaussianity is likely to be significant 
only for relatively small angular scales. 

It should be emphasized that the mechanism described previously is deeply rooted in the quantum-mechanical nature 
of the gravitational field. The observable quantities calculated in this framework arc always proportional to the Planck 
length. In other words, if observations confirm the full set of inflationary predictions then the fact that ST/T ^ would 
be a direct observational consequence of quantum gravity. 

The quantum-mechanical origin of the anisotropies in the framework of inflation raises also profound problems of 
interpretation. One should not think that these problems are purely theoretical. On the contrary, they have consequences 
with regards to the experimental strategy that one should follow in order to extract as much informations as possible 
from the data. The fluctuations in the CMB effective temperature are linked to the perturbed metric as shown for the 
first time by Sachs and Wolfe. Therefore, the fact that the perturbed metric is an operator implies that the primordial 
fluctuations in the temperature must also be considered as a quantum operator. The observables are often expressed as 
n-point correlation functions of the operator A(e) = 5T/T(e) in the arbitrary state \^>) 



e n (ei.---,e B ) = <*|A(ei)-"A(e n )|*>, 



(1) 



where e^'s are arbitrary directions on the celestial sphere. In the following, we will also use the notation £(ei-e2) = 
£2(61,62). According to the postulates of quantum mechanics, the previous theoretical predictions should be confronted 
to experiment in the following way. The same experiment should be performed N times giving each time different 
outcomes qi. If the quantity (1/-/V) X^j=i Qi E oes to the corresponding quantum expectation value when N goes to infinity, 
the theoretical prediction is said to be 'compatible with experiment'. This is the core of the problem in cosmology: we only 



have access to one realization, i.e. one map of the CMB sky and that means that N is fixed and equal to one. Therefore, 
the question arises as to how we can verify the theoretical predictions of the theory of quantum-mechanical cosmological 
perturbations. This is a way of stating the cosmic variance problem. It is a fundamental limitation in the sense that it 
remains even when other limitations like instrumental errors or low angular resolution have been fully mastered. 

The usual method to deal with this problem is to replace quantum averages with spatial averages over the celestial 
sphere. Suppose we wish to measure £ ra (ei, • • -e n ). (Of course, the discussion could also be applied to quantities other 
than correlation functions). The first step is to introduce a new operator, the estimator £(£ n ) of £ n , defined as 

= J ■ • ■ / dO x • • • dQ n £(£„)(ei • • ■ e„)A(ei) • ■ ■ A(e„), (2) 

where E(£ n )(ei ■ ■ -e„) is a weight function to be determined. Clearly, £(£„) is defined through a spatial average. The 
second step is to require that the estimator is unbiased, i.e. 

<tt|£(£„)|*>=*». (3) 

In general, this restricts the class of functions E(£ n ) allowed. The fact that the mean value of the estimator be equal to 
the quantity we are seeking does not guarantee that each outcome will be for sure . The third and final step is then to 
find the function E(£ n ) such that the variance (squared) of £(^ n ), i.e. 



T 2 



°i«„) = (^)£(Cn)) - (£ (£„)> 2 (4) 

be as small as possible, taking into account the constraint given by Eq. (^). The corresponding estimator is then called 
the best unbiased estimator. Mathematically, this requirement is expressed through the following variation equation 



<*|£(£n)|*>-£, 



= 0, (5) 



where one has introduced a Lagrange multiplier A which can then be determined from the previous equation and the 
constraint itself. Once we have A, we plug it into Eq. (^) and this completely fixes E(£ n ) an d hence the corresponding 
best estimator. In turn, its variance can now be calculated. If this one vanishes then we are sure that each outcome is 
and from one realization we can determine the n-point correlation function. In this case, A(e) is said to be ergodic, 
i.e. ensemble or quantum averages coincide with spatial averages. Unfortunately, one can show that this cannot be the 



case on the two-dimensional sphere (Grishchuk & Martin 1997). Otherwise, we have found the weight function E{^ n ) 
which leads to the smallest non-vanishing variance. If the variance is small enough, each outcome will be concentrated 
around the mean value and with just one realization we have good chances to get a reasonable estimate of the correlation 
function. The typical error made in considering that one given outcome is equal to the mean value is characterized by 
the variance of the estimator, see Fig. 1. 

All the above analysis performed for the best (quantum) estimators can equally well be reproduced in the case where 
the anisotropics are due to an underlying stochastic process although the former is generally physically best motivated. 
In that case, quantum averages (^| • • • |$) are just replaced with stochastic averages (• • •}. In the following, we will drop 
out the 'hat' symbol and consider that the different quantities are either operators or stochastic processes. In the same 
manner, we will denote an ensemble average by the symbol (• • ■) having in mind that this means either quantum or 
classical averages. 

Let us now describe the relevant quantities to estimate. It is convenient to expand the temperature fluctuations over 
the basis of spherical harmonics according to 



A(e)=5>ri7». (6) 



Cm 



5 





Fig. 1. — Sketch comparing the variances of an arbitrary estimator of £ n with the best estimator of the same quantity. The widest 
distribution does not permit an accurate determination of § n due to its large variance whilst the narrow distribution corresponds to the best 
(unbiased) estimator and possesses the smallest possible variance. As a consequence, one given realization (i.e. our CMB sky) will most 
probably be closest to the mean value in the second case rather than in the first one. 



This equation assumes a complete sky coverage. Implementation of the method for the incomplete (galaxy-cut) sky can 
be performed by using the basis introduced in ( |G6rski 1994 ). Once a specific model is given, the statistical properties of 
the a™'s are determined. Since A(e) is real, the a™'s must satisfy aj 1 * = (— l) m aj m . Without restricting the generality 
of the underlying physics, the first three moments can be written as 



/ tn\ n / Tiii rri2*\ n X £ I ~vn\ n 

{a e ) =0, {a^a^ ) = C^d^^ 2 d mim2 , [a ti a 



m 2 m 3 \ _ fix 
a e 3 I ~~ Ul 



«2 t3 \f 
rri2 m 3 ' ^ 



(7) 



where ^ 2 is a Wigner 3j-symbol. The second equation ensures the isotropy of the CMB. The quantity (a^a^ 2 *) 
is the second moment of the a™'s and Ct is usually called the angular spectrum. In the third equation, the quantity 
(a^ 1 a^ 2 a^ 3 ) is the third moment while Ci 1 g 2 £ 3 is called the angular bispectrum. For l\ — £ 2 = ^3 = A this quantity is 
generally written as Bi = Cm. The presence of the Wigner 3j-symbol guarantees that the third moment vanishes unless 
nil + fn% + 777,3 = and \i{ — £j\ < Ik < &% + £j. Moreover, invariance under spatial inversions of £3 implies an additional 



'selection rule' ( |Luo 1994 |Gangui fc Martin 200C| ), h + £ 2 + h = even, in order for the third moment not to vanish. 
Finally, from this last relation and using standard properties of the 3j-symbols, it follows that the angular bispectrum is 
left unchanged under any arbitrary permutation of the indices 

We will need the higher moments as well. Since departures from Gaussianity are expected to be small (specially on large 
angular scales), higher moments will be calculated in the mildly non-Gaussian approximation. Within this approximation 
we can write a™ = aj 1 ^ + taj 1 ' 1 ' + 0(e 2 ) where a™*- ^ is a Gaussian random variable and the expansion parameter e 
is small. In the following, each moment will be calculated to the first non-vanishing order in e. For example, the fourth 
moment yields (a™ 1 a^ 2 a™ 3 a ^ 4 ) = (a™ 1 ^a^ 2 ^a™ 3 ^a^™ 4< "°' ) ) + 0(e). As a consequence, the connected fourth moment 
can be neglected because it is of higher order than the Gaussian part. The 'w' label will be dropped out hereafter. 
Therefore, in the mildly non-Gaussian approximation we can write 



— {~^) m2+rnA Ci 1 Ci 3 5t 1 t 2 5 mi - m2 5t 3 t i 5 ms - mi + Ci 1 Ci 2 8g 1 (, 3 5 mim3 5t 2 t i 5 m2mi 

77117714 ^t%tz ^7712777,3 • 



(8) 



The fifth moment could be determined in a similar way but we will not need this quantity in the following. Finally, the 
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sixth moment can be expressed as 



(a™ 1 a™/ a™/ a™ 4 * a™ 5 * a™ 6 *) « (a^a£ a )(a£ 8 a^*)(a£ 5 *a£"*) + 14 additional permutations . (9) 

Although the explicit expression is not particularly illuminating, the last equation will be useful for the calculation of the 
variance when dealing with the three-point correlators below. In particular, one can write 

/ mi ni2 nt3 mi* ni2* m3*\ n n n i o/^3 c /r c \/r i X "\ 

\ a £i a £ 2 £j °^l a f-2 a i'i I ~ ^(1^(2^(3 + Z ^l 1 "eie 2 i3\ (, rn 1 m3 + 0m 1 -m 3 )(<Jm 1 m 2 + Om 1 -m 2 ) 

+ Ci> 1 C'f 2 8i 2 i 3 {8 m2m3 +<5m2-m 3 ) "I" ^4^4^Wl (^miras +^mi-mj) 

), (io) 

where the symbol (^^^ vanishes unless £\ ~ £ 2 — £3 in which case it is one. This equation coincides with Eq. (24) 



of ( Luo 1994 ) provided the undefined symbol S mim2m3 fi written in that work has the meaning i mim2m3i o = (S mim3 + 

&m\ — 1713 ) (^mim^ ~1~ ^mi — 7712 ) • 

As we mentioned in the Introduction, in this article we are mainly interested in finding the best unbiased estimators for 
the two following quantities: the third moment (a™ 1 a^ 2 o™ 3 ) and the angular bispectrum Cg 1 t 2 e 3 . These are related to 
the three-point correlation function £3. It is clear that, as mentioned above, the very same method could also be utilized 
for the computation of with n arbitrary. Before addressing this question, however, we will first treat the analogous 
quantities related to the two-point correlation function, namely <^a™ 1 a^ 2 *^ and Ci, the main purpose being to illustrate 
concretely the tactics presented above in a case where everything can be calculated easily. This will be used as a guideline 
for the case of the three-point correlators. 

3. TWO-POINT CORRELATORS 
3.1. Best estimator for the angular spectrum Ci 
The definition of the estimator £(Ci) is given by our general prescription, see Eq. (U) 

£(C e ) = J J dn 1 dn 2 E i (e 1 ,e 2 )A(e 1 )A(e 2 ), (11) 

where E (ei,e 2 ) is the weight function. The angular spectrum Ci is a real quantity and its estimator £(Ci) must also be 
real. Therefore, the weight function can be taken real. From the previous definition, it is clear that the antisymmetric 
part of E e (ei,e 2 ) does not contribute to the estimator. Then, we can replace i? (ei,e2) in £(Ci) by its symmetrized 
expression i5|(ei, e2) = (l/2)[E e (ei,e 2 ) + E e (e 2 ,ei)]. At this stage, two methods can be applied. Either we work directly 
with the weight function or we expand it over the spherical harmonics basis and try to determine the coefficients of this 
expansion. Clearly, both paths are equivalent and can be followed for any estimator. Here we employ the second method, 
leaving the first one for the determination of the second-moment estimator considered in the next subsection. Therefore, 
we write the weight function as 

£|( ei ,e 2 ) d f\ e 2 ^{ ei )Y^{e 2 ). (12) 

£imi e 2 m 2 \m 1 m 2 J 

The reality and symmetry properties of the weight function Eg imply that the complex coefficient of the expansion must 
satisfy, respectively 

df* s = (-l) mi+ma d/ -, , df , =df , . (13) 

{L\ll} { {XX} {11} 



Inserting the expression of the weight function (12) into the general definition of the estimator ([Ll]) and using standard 
properties of the spherical harmonics, one gets 



(14) 

Our first move is now to require that (£(Ce)) = Cg. Using the second of Eqns. ([?]), we find that the coefficients d must 
fulfill the following constraints 

" (15) 



mi \ _ mj mi J 



All the estimators satisfying this condition are unbiased. However, it is clear that this does not completely determines the 
estimator but just a class of estimators. Our second move is to calculate the variance. Using Eq. (|J), a straightforward 
computation gives 

^(C)= 2 E E d r\ «,l d ( (l i2 \ c ^- ( 16 ) 

This quantity is obviously positive. From this expression, one sees that the imaginary part of the coefficients d only 
increases the variance. Since a vanishing 3(d) satisfies the constraint equation, we can consider that the cfs are real. 

Our third move is to minimize the variance taking into account the constraint. For this purpose we introduce a set of 
Lagrange multipliers A| (i.e. one Lagrange multiplier per constraint since £ must be seen as a fixed index) and require 
that 



^(o+£*£(£(- 1 ) mid { „ o 



0. 



(17) 



tl v mi 

The definition of the variation 8 must respect the symmetry properties of the coefficients d\ we take 
5d' 



1 f 1 
— ( 

2 

1 



w mim 1 u m2m 2 ~ \ ± ) u m\ — m^ u m 2 — m 2 



8 ,8 , + (-i) m i+ m 28 1 8 



(18) 



Although it is not compulsory to use this equation, since a naive definition of the variation would lead to the same final 



result (Grishchuk & Martin 1997), it is nevertheless interesting to utilize it as a warm up for what will be done for the 
angular bispectrum. The variation leads to the following relation between the coefficients d and the Lagrange multipliers 

4dr € . . >C<A + (-l) roi A| 1 ^ 1 - m ,=0. (19) 



We see that the choice of the A s is not free; it is fixed by the variation itself. Using the constraint in this equation, 
we find \\ = —A\Cf/(2lx + l)]df . Having determined what the Lagrange multipliers are, the problem is completely 
solved. It is now sufficient to re-introduce this value for Af in Eq. (|l9|), get the coefficients d and, from this, also the 
best unbiased estimator: the weight function is then given by E$ Bcst( e i> e 2) = (l/47r)P^(ei • e2) and the estimator itself 
by (see also |Grishchuk fc Martin 1997| ) 



2f 



1 1 
— £ a 



m m* 
I a £ ■ 



The variance of this estimator is the well-known 'cosmic variance' 

2 _ 2Cg 



(20) 



(21) 



One remark is in order here. The cosmic variance is usually obtained in the following way: the previous estimator appears 
naturally from Eqns. (0) and it is usually assumed that the a™'s are Gaussian random variables. In this case, the estimator 
( |2C| ) has a x 2 probability density function and from this the cosmic variance can be easily recovered. The proof presented 
above is by no means equivalent to this naive derivation. There are many unbiased estimators and, a priori, nothing 
guarantees that the simplest one is the best one, i.e. the naive derivation is not sufficient to prove that the estimator 
( |20| ) is the best one. This can be proven only along the lines described above. Moreover, we do not need to assume that 
the a™'s obey a Gaussian statistics. Only the mildly non-Gaussian assumption is necessary for the calculation of the 
four-point correlators. 

3.2. Best estimator for the second moment C a 

In this section, we discuss the best estimator for C a = (a^ 1 a^ 2 *') . Greek letters will always be employed for collective- 
index notation, like a = ^ \ or a' = ^, \. Unlike in the foregoing subsection, we find the weight function 
directly without expanding it over the spherical harmonics basis. This method is closer to the one used by Heavens 
(1998). We also show that the best estimator of C a cannot be deduced from the best estimator of the angular spectrum 
Cg obtained in the previous subsection as one could naively think. 

Let us start with a couple of definitions. First, the quantity 

^(ex,e 2 )=y^(ex)y^*(e 2 ), (22) 

which is complex and non-symmetric with regards to the position of the complex conjugate symbol '*'. This last property 
comes from the definition of ^a™ 1 a^ 2 *) itself and it would be the case for any even-point correlator. Then, regarding 
complex conjugation, there is a slight difference between the even- and the odd-point correlators. The real part of 
i? a (ei,e 2 ) will be noted 



i?£(ei,e 2 ) = ~ 



Y e ^(e 1 )Y^*(e 2 ) + Y^*(e 1 )Y e ^(e 2 ) 



(23) 



-Rj£(ei,e2) is not symmetric under a permutation of the two directions. It satisfies i?^(e2,ei) = R^(ei,e 2 ), where the 
index a is defined by a = { ^ ^ }. The symmetries in the indices and in the directions will play an important role in 
what follows. Using the previous definitions, we can introduce a quantity which is symmetric both under a permutation 
of the two directions and under a permutation of the columns of internal indices in a 



i^(ei,e 2 ) = i 



Y™i( ei )Y™**(e 2 ) + r™ 1 (e 2 )r; 2 " 2 *(e 1 ) + Y^*( ei )Y^ (e 2 ) + Y^*(e 2 )Y^ (e x ) 



(24) 



i.e. we have i?g(ei,e 2 ) = i?g(e 2 ,ei) = i?g(ei,e 2 ) = J?g(e 2) ei). Finally, we will also employ a symmetrized Kronecker 
symbol defined according to 

^ a = 2 {^hi' 1 ^t-2t'^m 1 m'^m2mi 2 + ^i^^sf^mim^mamiji (25) 

which is left unchanged under a permutation of the indices a and a' and also under permutations of the columns of each 
collective index separately. 

Let us now turn to the computation of the best unbiased estimator; the general definition of an estimator can be 
expressed as 

£(C a )= J J dn 1 dn 2 E a (e 1 ,e 2 )A(e 1 )A(e 2 ). (26) 

The quantity C a is unchanged if we permute the columns of indices in a, C a — C & with a given above. Being an estimator 
for C a it is then natural to assume that £(C a ) possesses the same property, £{C a ) — £{C & )- Looking at the definition of 



-E a (ei,e 2 ), Eq. (pq), we easily see that the weight function will also satisfy -E a (ei,e2) = E a (ei,e2). Moreover, we take 
E a (ei,e2) symmetric under a permutation in the directions ei and e 2 , i.e. Eg(ei,e2) = Eg (e2,ei). 

Now, we require that the estimator £{C a ) be unbiased, which implies that the following relation must be fulfilled 



dn 1 dn 2 E^{e 1 ,e 2 )R a (e 1; e 2 ) =6%° 



(27) 



We have required the presence of Sg a in the right hand side of the previous equation to respect the symmetries in a of 
the weight function. In the previous equation the Eg and R a are a priori complex. However, one can show that it is 
possible to work only with a real weight function, 25g*(ei,e2) = J5g(ei,e2) and with the R^(e±,e2) defined above: the 
imaginary contributions would just increase the variance. (One could have also chosen to work with a pure imaginary 
weight function and the imaginary part of R a ). Therefore, the constraint can be written as 



d0 1 dfi 2 ^(e 1 ,e 2 )i?g, (ei,e 2 ) 



(28) 



where -Eg has been taken real. Let us now calculate the variance of the estimator £(C a ). Using Eqns. (g), we easily find 
that 

(29) 



2 

a £(C a 



2 / J J J dQ 1 dQ 2 dn 3 dil i Ei(e 1 ,e 2 )Ei(e a ,e i )C(e 1 -e 3 )^e 2 -e i ) 



Our next step is to minimize this variance under the constraint given in Eq. (|2q). For this purpose, we introduce a set of 
Lagrange multipliers A" and require that 



dn 1 dn 2 ^(e l! e 2 )^(ei,e 2 )-C 



= 0. 



(30) 



At this point the precise meaning of the variation symbol 5 matters. Before performing the variation, let us recall that 
the symmetries of the weight function must be respected; hence we will have 



SEgjeuej) = 1 
5El(e k ,e e ) 2 



5(ei-e k - I) 5(e r e e - 1) S(e i -e i - 1) £(e,-e fe - 1) 



2tt 



2tt 



2tt 



2tt 



°s ■ 



(31) 



The 27r's in the denominators come from the fact that, while the direction is expressed in terms of the corresponding 
spherical angles like = (9i,tpi) and then d£li = dcos8id(pi, there is an extra 2ir factor in 5(ej-efc — 1) = 2ir8(cos6i — 
cos9k)S((pi — ifk)- As a result of the variation, we obtain the following equation 



dn 1 dn2^(e 1 ,e 2 )C(ei-e 3 )^(e2-e4) 



-E 

4 ^ 



Kx'Rs (e 3 ,e 4 ) 



(32) 



The quantity Rg (e3, e 4 ) appears naturally as a result of Eq. (|3lj). The previous equation should be compared with Eq. 
( JTg| ) of the previous subsection. The result of the variation is a relation between the weight function and the Lagrange 
multiplier. Our aim now is to get an explicit expression for the weight function. This can be done by using the inverse 
two-point correlation function £ _1 which satisfies (see also the Appendix) 



1 < !<>,;"•: (\-<!, ii 1 {e r e k ) = 5(e r e fc - 1) 



(33) 



In the case of the three-point correlator, this definition leads to subtleties which will be examined in detail in the next 
section. Multiplying Eq. ( |32| ) by £ _1 (ei' • es)^ 1 ^' • 64), integrating over directions e3 and e 4 and relabelling indices, 
we arrive at 



£s(ei,e 2 ) 



1 1 



4 (2?r) 



2 L^i 



J yd0 3 dfi 4 A^r 1 (ei-e 3 )r 1 (e 2 -e 4 )i?g'(e3,e 4 ) 



(34) 



This equation for the second moment is the analogous of Eq. (21) of (Heavens 199S ) obtained for the third moment. It 
is clear from the previous section that, at this stage, our final goal has not yet been reached. The weight function is still 
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expressed in terms of the Lagrange multipliers. The correct way to proceed is to remove the latter using the constraint 



given by Eq. (|28|) as it was done in the previous subsection. Then, we first multiply Eq. (34) by the quantity Rg (ei, e2) 
and, after having integrated over solid angles fix and fl 2 , this leads to the equation that the Lagrange multipliers must 
satisfy 



(35) 



This equation is the analogue of the equation Af = —4[C^ i /(2£i + 1)] <5j of the previous subsection. Here, the difference 
is that a combination of Lagrange multipliers with different indices appears rather than the Lagrange multiplier itself. 
However, we can reconstruct exactly this combination in the right hand side of Eq. (|34|). Indeed, we just have to multiply 
each side of Eq. ( |35| ) by Rg (ei,e2) and perform the sum over a". Using the symmetry properties of Rg (ei,e2) we 
obtain 

4C^ 2 i^(ex,e 2 ) = ^A«„i^"( ei ,e 2 ). (36) 



We now insert this equation in the right hand side of Eq. (34) to remove the Lagrange multipliers and find 

£ S Q Bc S t(ei,e 2 ) = ^2 C ^J J dn 3 dn 4 C 1 (e 1 -e 3 )C 1 (e 2 -e i )R^e 3 ,e 4 ). (37) 

To go further, we express in Eq. (p7|) explicitly in terms of spherical harmonics (see the Appendix) which yields 
Bost( e i' e 2) — ^s( ei ' e2 )- Now, we just plug this result into Eq. ( ]26| ) to get the final explicit expression of the best 
unbiased estimator £-Bcst(C a ) 

c In \ ( mi* m? I mi mi* i /oo\ 

tBest(Ca) = -I a^ 1 a l2 2 +0^0^ J . (38) 

Some remarks are in order here. First, despite appearances, one cannot deduce £ , Bcst(C Q ) from fBcst(C£ 1 ). It is clear that 
the following equation holds 

(£(C a )) = (£(Ce 1 ))5e 1 e 2 5 mim2 . (39) 

However, from this equation we are not allowed to conclude that £Bcst(Ca) = £Bcst(Ci 1 )5e 1 e 2 5 mim2 . Therefore, knowing 
one of the best unbiased estimators does not allow us to infer the other best one. To be specific, if one assumed the 
previous wrong relation, from Eq. ( pb] ) one would get 

WC«) = 2JTi( «r<)<Wmim 2 (false), (40) 

and although this estimator is unbiased, it is not the best one. The second remark is that although the estimator given by 
Eq. (|38|) could be naively regarded as a trivial one, this is not the case. Indeed, the estimator of Eq. (^) is as trivial as 
the actual best one. The moral is then that there exist simple choices which lead to the wrong answer. The only reliable 
method in the problem of minimizing the variance is therefore the one exposed above. 

4. THREE-POINT CORRELATORS 

4.1. Best estimator for the angular bispectrum C^ 2 ^ 3 

In this section, our aim is to determine the best unbiased estimator for the angular bispectrum C^e^ defined in the 
third of Eqns. ([?]). According to our general prescription, the most general definition reads 

£{C llM3 ) = J J J dn 1 dn 2 dn 3 E e s Ae >(e 1 ,e 2 ,e 3 )A(e 1 )A(e 2 )A(e 3 ). (41) 
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As in the case of the two-point correlators, the weight function also possesses the properties of being real and symmetric 
under arbitrary permutations of directions e^. In addition, like Ci x t 2 i 3 ^ the weight function satisfies 2 3 = E s 2 as 
well as for any other arbitrary permutation of the indices 1%. We follow similar steps as for the angular spectrum and 
therefore we choose to expand the weight function on the basis of the spherical harmonics. Then, as in Eq. (|l2|), we write 

£^ 3 (ei,e 2 ,e 3 ) = ^ E E d f \ \ \ i (e^ (e 2 )Y^ (e 3 ) . (42) 



^ f-; 7 4 4 4 \ 4 



The properties of the weight function imply that the coefficients d must satisfy equations similar to those given in Eqns. 

L m l m 2 m 3 J I ~ m l _TT1 2 ~ m 3 J L m l m 2 m 3 J I m 2 m l m 3 J 

where the last relation is in fact valid for arbitrary permutations of any two columns of the collective subindex. Like the 
weight function, d is also left invariant under arbitrary permutations of indices l\ (not primed). The estimator can be 
expressed in terms of the coefficients d and the a™'s only: inserting the expansion of the weight function in the above 
expression for the estimator and using standard properties of the spherical harmonics one obtains 

£(cw 3 ) = e e e *;\ \ \ r-t^si- ^ 

i\m' x t' 2 m' 2 t' 3 m' 3 \ m ; m ' 2 m > 3 j 

In practice, CMB observational settings are devised such that both the monopole and the dipole are subtracted from 
the anisotropy maps. This means that the coefficients d in the last equation are only non-vanishing for indices £' i > 2 in 
the collective subindex. Moreover, the coefficients d satisfy ^1+^2+^3 = even. We must now require that our general 
estimator given by Eq. ( p4| ) be unbiased, i.e. (£(Ci 1 e 2 e 3 )) = Cg t g 2 g 3 . This forces the coefficients d to fulfill the following 
constraint 

,„' ,„'. ) ,„'■ ,„" ,-' ( \ ' n ' 2 m 3 

1 2 3 L l 2 3 J V 

where we have defined a new symmetrized Kronecker symbol, this time for the £ multipole indices only, as follows 



tit 1 - 1 / \ 

S s 3 = — I Si 2 e' 2 5i s i> s + 5 additional permutations J. (46) 

It is easy to check that the constraint equation satisfies the conditions imposed by Eqns. ( pl3| ) on the coefficients d. In 

hi' 

particular, let us justify the presence of the symbol S s 3 . Using the previous properties for d, relabelling the indices 
m\ «-> 77?2 in Eq. fllq ) and finally noting that £[ + £' 2 + £' 3 = l\ + £ 2 + ^3 = even, which allows us to permute any two 
columns of the Wigner 3j-symbol, one verifies that the left hand side of the constraint is invariant under £[ <-> £' 2 . The 

til'- . a — I 

same applies for any pair of £ multipole indices and this explains the presence of the symmetrized S s 3 in Eq. ([45| ) . We 
see from this that all coefficients d that do not satisfy £[ + £' 2 + £ 3 = even do not enter the constraint. We will show below 
that these terms only increase the variance and as a consequence one can take them equal to zero. In particular, note 
that this is the case for a coefficient d with £[ = £' 2 and £' 3 = odd. This property will turn out to be useful in what follows. 

We are now in a position to calculate the variance of the estimator. Looking at Eq. we see that this requires the 
computation of the sixth moment of the et™'s, see Eq. (g). After having made use of the properties of the coefficients d 
and rearranging the resulting 15 terms into two groups, straightforward algebra yields 



([£ (CW 3 )f ) = E E E C 4 C 4 C * 



6d n «2 (a. d ti *2 ^ Q(_ 1) m' 1 +m' 2d «i ^ e s * j H '2 e s 

j l \ f '2 e 3 \ j e 'l £ 2 >3 X ( S e '\ S X { ^ e 2 r e 3 X 



(47) 



12 



The square of the variance of £(Ci t i 2 i 3 ) is given by 



(48) 



Following the discussion in §2, the term (\£{Ct l t 2 t 3 )f') is of order e° whereas the lowest non- vanishing order of (£(£^^3)) 
is e 2 . Therefore, the latter one will not enter the minimization procedure and the variance squared will be written as 



£(Ci 1 i 2 e 3 ) 



( x [£(Ci 1 e 2 g 3 )) 2 '). However, let us notice that this does not occur in the case of the two-point correlator. Indeed, 
as we have seen, in that case both terms contributing to the square of the variance are of the same order in e. Then, it 
follows that Eq. ( fi"7| ) corresponds to Eq. ( |l6| ) in the last section. 

Let us now examine the structure of the variance in more detail. The term 6c?* g? oc 5ft 2 (c?) + 3 2 (g?) in Eq. p7| ) is 
analogous to the one in Eq. jl^). However, here there is another contribution, the 9c?* d term, which will play a crucial 
role in what follows. The imaginary part of this term of course vanishes, as the variance must be real. However, the real 
part of it contains a contribution like 



it in'-. 



tL in' 



PL mL 



(49) 

where the quantity [. . .] 2 depends on the indices ^1,^2,^3 and i^m'^ and is strictly positive or zero. Of course, there is 
a similar term coming from the real part of d. Thus, we see that the various contributions of the imaginary part of the 
coefficients d to the two terms, 6c?* d and 9g?*c?, only increase the variance. Since we know that a vanishing imaginary part 
does satisfy the constraint Eq. ([iB]), it can be disregarded in the sequel. Therefore, Eq. (|47]) can then be written solely 
in terms of real coefficients d as follows 

2 



£{Ce l e 2 t 3 ) 



E E E c <i c 4 c <i 



6[d fl 12 ' a +9(-l)< +m '^d ei ' 2 ' 3 d' 1 ' 2 * a 

L m l m 2 m 3 J L m l ~ m i m 3 J I m 2 ~ m 2 m 3 J 



(50) 



Our next move now is to minimize this variance with respect to the coefficients d, taking into account the constraint of 
Eq. © 

£ 1 1' 2 1' 3 
m' x m' 2 m' 3 

This equation is the analogous to Eq. (|l7|). As before, one needs to give a concrete meaning to the symbol 6 of the 
variation. Its definition must respect the symmetries of the coefficients d and hence we take 



' 3 \ m [ m> 2 m ' 3 J 



0. 



(51) 



Sd 



( S % % ) 

L m l m 2 m 3 J 



Sd 



t-l «2 «3 



( %, % ''I ) 
I m l m 2 m 3 J 

Then, Eq. @ leads to 



— — S ~3f P" $P' P" $f f" 



5 terms L (52) 



/ /' /" /" 

UCeCeCed 1 } 1 ? 1 *,, +X i //,f 1 2 3 

l l l 2 l 3 \l',t'„l'„\ \ I I I 

1 m\ m 2 m 3 



1*2 '■?, 



{&9 



ii-i^^Eui-irt: , ,, j 



6(-l) m 3 C ,,C ri ^^_ m ^C £ (-ird flf ;^ +6(-ir'xC, i C^^5 m ,_ mi ^C,(-ird^;^ =0.(53) 



(')» 



This formula, together with Eq. ([15[), form a set of equations which completely determines the best unbiased estimator. 
We see the complicated structure of these equations. The last three terms come from the 9c?c? term in the variance and 
are not present in the case of the angular spectrum. 
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From this last equation and using the constraint Eq. ( |45| ) we can get the general expression for the Lagrange multipliers. 

Thus, we multiply Eq. ( |53| ) by the appropriate 3j-symbol and we sum over the three indices m^. The first term is exactly 

lit'- 

the constraint and produces a S s 3 . Using the fact that a triple sum over the mi's of the squared of a 3j-symbol gives 
unity, the second term yields the Lagrange multipliers themselves. Finally, the last three terms vanish: indeed, after 
straightforward manipulations one generates a term like 



\m -m J "h 



£■> £? 



■I 

-771 



(54) 



for the third term of Eq. (|53j) and analogously for the last two ones. As we mentioned previously, the coefficient d in Eq. 
d54| ) vanishes unless £' 3 is even. In this case, recalling the identity (Mollerach et al. 1995) 



(-1)' 



£[ £[ 







= (-i) £ ' lv ^TTfco, f a 



even, 



(55) 



we see that there will only be a non- vanishing term if £' 3 — 0. But, the corresponding d is zero because £' 3 < 2 and therefore 
terms of this kind do not contribute to the Lagrange multipliers. Then, these are given by 

A«»?»? = — YICojCii Cei st 3 . 



(56) 



Plugging this into Eq. (p3|) , one has 



12Cf'Cf'Cf 



iA [ 4 £ 2 ^ 



! e 'i e 2 e 'A 
h'i n 'i"'J 

+ Q{-l) m '*C l ,C l , l 5 V > 3 5 m>2 _ 



l 3 X] 

lm 



m\ m, mi 



c e (-i) m d e ; e2e3 



+ et-ir'^c^^ 



rnU—m 



{m-mm'J lm 

This is the final equation to be solved in order to determine the best unbiased estimator. A solution is 



E^rdp ^|=0.(57) 



This leads to 



j «1 ^2 '3 

{ l \ s s \ 

I m l m 2 m 3 J 



£Bcst{Cl 1 l 2 l 3 ) — ^2 



&-t $o 



'-2 

m' 2 



3 

m', 



*3 

mi 



m-,^ m 2 m 3 



(58) 



(59) 



This is the main result of this subsection. 

Given that we now know the best unbiased estimator for C^ 2 ^ 3 , one can compute its variance, the smallest one amongst 



all possible estimator variances. In ( Gangui fc Martin 2000 ) we have already calculated it and reads 

a SBest(c ei e 2 t 3 ) = C h c t 2 Ci 3 (^ + <W 2 + 5t 3 t 3 + <Wi + 2 S ei £ 2 Se 2 i 3 ). (60) 

In the same reference a plot of this variance for low order multipoles can also be found. This is what one could dub (the 
square of) the 'bispectrum cosmic variance' in perfect analogy with <t| b , C s = 2Cf/(2£ + 1), which is (the square of) the 
variance of the best unbiased estimator for the angular spectrum, commonly known as the 'cosmic variance'. 

Let us conclude this subsection by comparing our results with those recently appeared in the literature. An estimator 



restricted to the diagonal case £i = £2 = £3 has been proposed in (Ferreira et al. 1998, see also Magueijo 1999 for an 
extension of their analysis) for Bg = Cue and reads 

-3/2 




E 




(61) 
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In that work, the aim of the authors was not to seek the best estimator, but to use Eq. (61) to analyse the non-Gaussian 



features of the 4-yr COBE-DMR data. It is easy to see that their estimator does not satisfy the constraint (Eg), i.e. the 
estimator is biased. This is due to the presence of the overall prefactor in front of the triple sum in Eq. (|6l|). However, 
as we have proven above, getting rid of it produces the best unbiased estimator £&est(Ct 1 e 2 e 3 )) Eq. (|59|). 

4.2. Best estimator for the third moment B a 
We now seek an estimator for B a = (a™ 1 a^ 2 a^ 3 ) where, as in the last section, it is convenient to define a collective 



index a = { m 2 m 3 }. This question has already been addressed in ( Heavens 199S ). As we did with the second 
moment, our starting expression for an unbiased cubic estimator £(B a ) of B a will be in the form 

£(B a )=J J J dn 1 dn 2 dn 3 E a (e 1 ,e 2 ,e 3 )A(e 1 )A(e 2 )A(e 3 ) (62) 

The goal is to find the weight function ^"(ei, e 2 , e 3 ) that minimizes the variance of the estimator. As above, the quantity 
B a is unchanged if we permute arbitrary columns of indices in a: B a — B a , where for instance a = { ^ ^ ^ }. £(B a ) 
has the same properties as B a and then it follows that £{B a ) — £(B a ) for any column-permutated a. This implies that 
the £ Q (ei,e2,e 3 ) satisfies E a (ei, e 2 , e 3 ) = E a (e 1 , e 2 , e 3 ). 

Now, is E a (ei, e 2 , e 3 ) also symmetric under permutations in the directions ei,e2,e3 ? From its definition we cannot 
know, for these directions are integrated over in the above defining equation. Unlike the case for the second moment 
discussed before, here E a cannot be decomposed into a symmetric and antisymmetric parts. However, we can always 
write E a = Eg + something, and show that this last contribution to Eq. ( p2| ) vanishes. Therefore, there is no loss of 
generality in working with Eg (ex, e 2 , e 3 ) = ^[E a (ei, e 2 , e 3 ) + 5 terms ] which is symmetric under arbitrary permutations 
of directions e^. 

Demanding the estimator £(B a ) to be unbiased, (£(£!")) = B a , yields the first constraint equation that the weight 
function Eg (ei, e 2 , e 3 ) must satisfy 

'dr! 1 dr! 2 dr!3i?s Q (ei,e2,e 3 )i?^(e 1 ,e 2 ,e 3 ) = (5r', (63) 

where R a (ei,e 2 ,e 3 ) = Y™ 1 {e\)Y™ 2 (e 2 )Y™ 3 (e 3 ) and i?g (ei,e 2 ,e 3 ) is its real part. The form of R a comes from the 
expression of £3, viz. (A(ei)A(e 2 )A(e 3 )) = J2 a B a R a (ei, e 2 , e 3 ). The symmetrized Kronecker symbol can be written as 
Sg a = ^{Sg^i 6 mim i Sg 2 (i S m2m i Sg 3 £i S m3 m' a + 5 terms) as required to comply with the symmetry under permutations in 
the columns of a in Eg (e l5 e 2 , e 3 ). In the above equation, (e!,e 2 ,e 3 ) is clearly non-symmetric under a permutation 
of directions ei,e 2 ,e3. However, as above, we can define a symmetrized combination Rg(ei,e 2 ,e 3 ) = ^[R^(ei,e 2 ,e 3 ) + 
5 terms] (12 terms). Symmetrizing either in directions or in the columns of a in R^(ei 1 e 2 ,e 3 ) yields exactly the same 
Rg(ei,e 2 ,e 3 ). In the last equation and in what follows the weight function Eg is real for reasons similar to the ones 
exposed around Eq. ( p8| ) in the last section. 
Proceeding as in §3 and using Eq. (^), one gets 

a s(B") = J d ^i • • •/ dQ 6 Eg (e 1 ,e 2 ,e 3 )Eg {e i ,e 5 ,e 6 )£_(e 1 -e 4 ) 6£(e 2 -e 5 )£(e 3 -e 6 ) + 9£(e 2 -e 3 )£(e 5 -e 6 ) , (64) 

where, utilizing the symmetry of the coefficients Eg (e, , e 3 , e^) under e-direction permutations, only two types of £ products 
remain: first type, six terms where all the three £'s mix directions of the first and second Eg's and, second type, nine 
terms where only one £ [in the above equation, £(ei-e4)] does it. 

To minimize the remaining variance under the constraint (|63|) we introduce a set of Lagrange multipliers A a and write 

-2 \ ^ 



AS/ / / / dfi 1 dfi 2 dr>3^(e 1 ,e 2 ,e 3 )^(e 1 ,e 2 ,e 3 )- ( 5r =0. (65) 
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As already noted for the second moment, the symmetries of the weight function must be respected in the variation; 
hence we have 

-S(e 1 -e l - 1) S(e 2 -ej - 1) S(e 3 -e k - 1) 



6Eg(ei,e 2 ,e 3 ) = 1 



2vr 



2tt 



2tt 



5 terms 



cap 

°s ■ 



(66) 



Now, we vary Eq. (p5[) and, after some algebra, we come up with 



dfi 4 dfi 5 dfi 6 i|f(e r e 4 ) 6£(e,--e 5 )£(e fe -e 6 ) + 9£(e,"efc)£(e 5 -e 6 ) +5 terms j Eg (e 4 , e 5 , e 6 ) = A",i?g' (e i7 e i7 e fc ), 

a' 

(67) 

expression which is symmetric in the arbitrary directions ej,ej, efc as it should (note the presence of Rg ). 

We aim at getting an explicit expression for Eg (e4, e§, e$). A glance at the previous equation shows that we need to 
multiply both sides of it by the inverse of the correlation function £ _1 defined in Appendix. Concretely, we multiply Eq. 
( |67| ) by £~ 1 (ei-e^)£ _1 (e_ 7 '-e 7 '/)£ _1 (efc-efc') and integrate over directions e^e^e^; we get 



6 Eg(ei>,ej>,ek') 

, n(2-) 2 



dft 5 dft 6 £(e 5 -e 6 ) £ (e</ •e j >)Eg(ek', e 5 , e 6 ) + £ (e^-efe^Eg (e^, e 5 , e 6 ) + £ (e^-e^-Eg (e,v,e 5 ,e 6 ) 



da i d%dfi fe r 1 (e. i -e. i Or 1 (e J -e J v)r 1 (e fc -e fc O^A^i?^'(e 4 ,e,,e fc ) 



(68) 



As before, the 27r factors in the left hand side come from operations like 27r/(e^) = J" df2fc<5(ej ■ efc — l)/(e&), for an 
arbitrary function /(e^). However, we don't have an explicit expression for £g yet. We see that expressions of the type 
£(e5 •e6)f _1 (- • -)-Eg (-, es, eg) are the ones that prevent us from isolating the weight function. To deal with this, it is 
convenient to construct the combination JJ d^dQ^Eg (_, es, eg)£(e5 -eg). To reach this goal, we multiply both sides of 
Eq. ( |6"g| ) by £(ej> -e^) and integrate over directions e^and ew ■ This operation produces a divergence in the left hand 
side of this equation in a form of a Dirac function '5(0)'. In the continuous case, all methods lead to this unavoidable 
problem and although it has already appeared in the literature ( Heavens 1998 ), it has never been treated so far. That this 
divergence is a mathematical artifact we can see from the fact that, in practice, we never deal with an ideal experiment: 
the problem is solved when we take into account the fact that each different experimental setting is limited by a finite 
angular resolution. This is usually quantified in terms of an ^-dependent window function Wi (the circularly symmetric 



pattern of the observation beam in £-space), although more involved scanning techniques are also employed (White & 



Sredn .cki 1995 , Knox 1999 ). Then, only a finite number of multipolcs will effectively contribute to the correlation function 
and, as a result, the above mentioned divergence is regularized; indeed, we have 



C(erej) = 



2£+ 1 

47T 



CeWfPiiei-ej), 



which, upon using the expression for £ 1 given in the Appendix, leads to the quantity 



Afe = / d<> ; -(veX 'u V e,; = 27rYJ WfY e m ( ei )Y e m *{e k ). 



In particular, the previous divergence '(5(0)' now becomes 



E 



2£+ 1 



Wf. 



(69) 



(70) 



(71) 



This is a more realistic and finite object to work with in the case where two directions on the microwave sky coincide for 
a given experience. Notice that for an ideal experimental setting in which the window function Wg — > 1, or equivalently 
the beam a — > in the case of a Gaussian profile, A blows up. Since the terms of the type 9££ in Eq. (|6^) are not 
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Fig. 2. — Values of Aik (as defined in the text) with fixed direction e; at the center of the plot and direction ej. scanning a square of 
side 90° . We show the result for the COBE-DMR window function. The full width at half maximum of the plot is roughly 10° , of the same 
order as the resolution of the COBE-DMR detector. We can then expect the plotted 'spike' to select just one experimental pixel on the sky 
map, making the relation Aih ~ <5(ei-ej, — 1) a good approximation. Note that higher resolution experiments will yield a more peaked curve, 
although one would expect the goodness of the approximation to remain roughly the same. 



present in the case of the two-point correlators, this problem did not appear there. From now on, strictly speaking, all 
expressions should incorporate the window function. However, in what follows and for computational convenience, we 
will keep Aik ~ 5(ei-ei~ — 1) for =/= e^, a good approximation as we can see from Fig. 2. 

Endowed now with the above regularization method, we present the term JJ dfl^dQeEg es, eg)£(e5 • eg) in the 
following form 

dfi 5 dtt 6 £; s Q (e4, e 5 , e 6 )£(e 5 -e 6 ) = ^ J2 X «' III ^d^d^C 1 (e 4 -e^ -1 (e? -ek>)R%' (e*, , e f , e k , ) (72) 

a' 

with T> = 247r 3 [6 + 9 ^ 2 ^ (2 x 2tt + 47T.4)] -1 where, as expected, the factor A appears explicitly. Now, we just replace the 
six terms with prefactor 9££ in the left hand side of Eq. ( |67| ) with Eq. (|72| ) resulting in 

dfi 4 dn 5 dfi 6 ^(e l •e 4 )^(e J ; -e 5 )f (e fc -e^E^ (e 4 , e 5 , e 6 ) (73) 

= ^E A 2'{^s'(e i! e i! e )t ) 

a' 

- I J dflfdClk'^iej'-ek') R$ (e;, e r , e k >)£(ej -e fc ) + R$ (ej,e 3 v,ejfe/)f(efc-ei) + iZ§ (e fc ,e J v,e fc /)^(e i -e J ) |. 

This contains just one appearance of E a ; to get the weight function explicitly, we only need to multiply the equation by 
£ _1 ( e j -e i)£ _1 ( e .r e 2)£ _1 ( e fc' e 3) and integrate over directions e^e^e^. Finally, we get the expression for Eg in terms of 
the Lagrange multipliers A" 

£ s Q ( ei ,e 2 ,e 3 ) - __i_^^A« { J J J dn i d%dn fe r 1 (erei)r 1 (e J --e 2 )r 1 (e fe -e 3 )^'(e i ,e J -,e fe ) (74) 
-V J J J d^dftjdftfci^^e^eOr 1 ^^ 

Reached this point, we have an explicit expression for Eg , but still dependent on the Lagrange multipliers. This equation is 
well defined as it contains the renormalized quantity A. Within a particular experiment with a given resolution, the value 
A takes depends on what one means by two coincident directions. For example, for the COBE-DMR window-function 
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specification (a Gaussian beam with dispersion a = 3°.2) this yields roughly A ~ 158.5, including the quadrupole. It is not 
difficult to extend this to other scanning techniques. The previous equation is the analogue of our Eq. (M) corresponding 



to the second moment and also to Eq. (21) of (Heavens 199? ). In that article, a similar analysis is done but for the 
discretized CMB sky. Remark that no divergence appears in his case. Indeed, all relevant quantities are finite when 
evaluated for two directions pointing towards the same pixel. Our prefactor T> corresponds to 3/(2 + 32V) in that paper, 
where N represents the number of pixels in the map. Clearly, N — > oo when the pixel size goes to zero, as well as 
A — > oo when the window function Wi — > 1. At this intermediate step our corresponding expressions need not coincide 
because both depend on the particular regularization scheme used (be it discretization or usage of a window function). 
Despite appearances, we will show below that the final expression for the best unbiased estimator does not depend on 
these schemes. This cannot be inferred from Eq. ([74]) because we still need to remove the Lagrange multiplier. Unlike 
what was done in (|Hcavcns 1998 ), we now proceed further and express the weight function explicitly. Hence, we multiply 
both sides of Eq. (Q) by i?S (ei,e2,e3), where a" = { \ % * 3 „ }, then integrate over directions ei, e2 and e3 and, 

n, mj m 2 m 3 

upon using the constraint equation (|6^), we get 

l262 a " = - „ „ , ( \(-l) m "+ m 2+ m 3 \°i „ „ +\ a f „ „ 1 +5 additional permutations in index a") 

GOi'lCi'lC?" ll J l l 2 e 3 I J l l l 2 l 3 LJ J 

I ~ m l ~ m 2 ~ m 3 J \. m l Tn 2 Tn 3 J 



-1)' 



im 1 1 m - m 



(-l) m 3A Q r „ , +A Q f „ , 



1 — ► 2 
2^3 
3 — » 1 



(75) 



Eq. ( |7^ ) is the final algebraic equation that the multipliers must satisfy. For fixed a, a natural way to proceed would 
be to get an explicit expression for A", . Another way to solve the problem goes on a line analogous to the case of the 
two-point correlators: we just need to identify the complicated combination of Lagrange multipliers in the right hand side 
of Eq. ( [75| ) with the one in the right hand side of Eq. ([73]) [or, equivalently, Eq. ([74])]. In order to do that, we now 
multiply both sides of Eq. ( |75| ) by (e$, e^, e fc ), perform the six sums over the indices in a" and we end up with 



R'i{e l ,e J ,e k )Ct l Ct 2 C i 



Rr (ei,ej,ek)Ce»Ce>iCg 



jraa 



-V V 

36 



[AJ 



LM Im 



(-1) 

C L 



M 



(mM-JWj \ —m — MM\ \ -MmM \ 1 M-m-M j 1 Af - Mm J \ 



+ A" f 



+ (-l) m A Q r 



+ A a f 



+ (-l) m A a r 



LLP. 
-MM-m 



(76) 



This is the equivalent of Eq. (|35| ) . The aim now is to show that the combination of Lagrange multipliers in the right hand 



side of the previous equation is precisely the one which appears in the right hand side of Eq. (73). In the latter, let us 
express both i?g and in the second term in the right hand side in terms of spherical harmonics. After some algebra 
we get 

dfi4dft 5 dft6£(ei-e4)£(ej-e 5 )£(efc-e 6 ).Eg (e 4 ,e 5 ,e 6 ) 

I . T) f—T\M 



12 



\ a a ,R% (e i)ei) e fc ) 



°EE H) 



36 



LM Im 



C L 



K/"(ei)£(ei-e fc ) + Y e m (e 3 )ae k -e t ) + ^(e^fe,.^) 



\x a f , + (-l) m A Q r s +\ a ( , +(-l) m A Q r , +\ a ( , +(-l) m A Q r , 

ImAf — MJ 1 -m-MMj "I -MmM j |M-m-MJ W-MmJ 1 -MM-mJ 



(77) 



which, as advertised, yields the same combination of Lagrange multipliers of Eq. ([76j). Then, putting the last two 
equations together, multiplying by £ _1 and integrating three times, we finally get the weight function associated to the 
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best unbiased estimator 

£ sVt( e " e J. e *) = ^s( e ^ e ^ e fc) 
which implies that the best unbiased estimator itself is given by 



£ Best (B<*) = -(a^a 



(78) 



(79) 



This is the final answer and it is a new result. Let us make a few remarks. Firstly, this does not depend on A which 
shows that Eq. ([79]) is independent of the regularization scheme used. Secondly, Luo (1994) used the following complex 
unbiased estimator: £(B a ) = a™ 1 a^ 2 a^ 3 , although he did not claim it to be the best one. Thirdly, as for the two-point 
correlators, one cannot use Eq. ( [79] ) in order to infer ^Best (C^ 2 ^ 2 ). Indeed, promoting the equation 

h £2 £3 
mi m,2 rri3 

valid for the mean values of the estimators, to the following equation 

/ h i* e 3 

mi m 2 m 3 

valid for the estimators themselves, is an unjustified step. If, nevertheless, we used this false relation we would get 



/ mi 



m 2 „m 3 \ _ 
a £ 3 I ~ 



£ 



Best 



--Best 



(false) , 



(80) 



(81) 



&eBt(B a ) = 



which cannot be cast into i 



h £2 £3 

mi !Ti2 TO3 



E 



£1 £2 



m x m 2 



1)1, 



in--, 



m.. 



(false), (82) 
So, like for the two-point correlators, we see that one cannot 



infer the £Bcst(£> Q ) from £-Bcst(Ce 1 e 2 e 3 ) and vice versa. 

Endowed now with the best unbiased estimator, one can compute its variance squared, the 'third-moment cosmic 
variance', which reads 

- 1 

a e Beat (B") - 2 

+ C^C* Si 2 i 3 2S mi o8 m2 - m3 



fa 



Cf„Cf 



1 + <5mi0^m 2 | 5m30 



ci 5, 



+ 5n 



(3mim 3 



C£ 2 C e3 5e 3tl 



2^77120^7713 — 7711 ~t~ ^7713! 



+ Cg 3 Cj l 8i 1 i 2 2S m3 o5 rni - m2 + 5 mim2 + & mi -m 2 |- (83) 

For example, from this we can now compute the cosmic variance for the third-moment estimator £-Q est (B^) with j3 = 
{ 1 1 _ 2 } where £ = even and m ^ 0. This particular case is often treated in the literature, see e.g. (Luo 1994, Heavens 
1998). We find 

^BcstlB 13 ) ~ whereas the variance of the estimator used in (Luo 1994) yields Cg^g^) — 2Cf. Another 



} where \im\ \ m j\ f° r anv hj- With this choice we get ct^ b t (g 7 ) =Cf/2 



example comes from taking 7 = { m ^ 
whereas Luo (1994) obtains — ^e- Note that in both examples the results differ by a 1/2 factor. This can be 

traced back to the form of the best estimator in Eq. ([79]). The variance computed in (Luo 1994) is consistent with 
his choice of the estimator. However, unlike what is stated in that paper, this variance does not deserve the name 
'cosmic' because, as we saw above, this estimator is not the best one. As expected, the variance of the best estimator 
is smaller than the variance computed in his article. Since we have now the correct expression for the cosmic variance, 
its numerical value should be re-estimated. To be specific, let us take (3 = { 2 2 _ 2 2 }. The cosmic variance is then 

1.3 x 1CT 15 where we used T = 2.7 K and 



°s Beet (Bf>) = Cl /2 = (4?r/5) 



^^ns-Ps/^O 3 



-PS = 18.7/xK (Runn & 



White 1997 ). Although this figure is close to Luo's result (1.4 x 10 -15 ) cited after equation (32) in (Heavens 1998), this 
does not imply that the two variances are not different by a factor 1/2, as we have just seen. This might probably be due 
to a difference in the quadrupole normalizations. 
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5. CONCLUSIONS 

Optimized analyses of CMB datasets involve the use of appropriate methods in order to reduce the various uncertainties. 
In particular, the theoretical error bars due to the cosmic variance can be minimized by working with the method of the 
best unbiased estimators. In this article, we have applied this technique for the study of CMB non-Gaussian features. 
These are often characterized by means of the third moment for the a™'s or by the angular bispectrum Ce 1 e 2 e 3 - We have 
found the best unbiased estimators in both cases. These are the quantities that should be used in future data analyses and 
would be important for upcoming megapixel experiments ( Tegmark 1997 , Borrill 1999 ) like MAP0 and Planck Surveyor^ 



In addition to this, we have displayed both the angular bispectrum and the third moment cosmic variances, the smallest 
possible uncertainties attached to the bispectrum and the third moment, which would be present in any ideal experiment 
when all other sources of noise have been removed. 

ACKNOWLEDGMENTS 

A.G is member of CONICET, Argentina. 

6. APPENDIX: THE INVERSE TWO-POINT CORRELATION FUNCTION 

In this Appendix, we derive the exact expression of the inverse of the two-point correlation function defined 
according to 

J .io-uvrX ';cv<- = 5(e v e k - 1). (84) 
As usual, one expands the two-point correlation function on the basis of the Legendre polynomials 

e(e,~e,) = ]T ^tic/flKere,-). (85) 

i 

In general, £ _1 can be expanded on the spherical harmonics basis as follows 

^(ere,) = £ ]T ^w^e^Y/ 1 '*^)- (86) 



Im I' m' 



Our aim now is to determine the coefficients bwmm'- Using Eqns. (85) and (|86|), the completeness relation for the 



Legendre polynomials together with the addition theorem of spherical harmonics in Eq. (j§J) , one gets 

E Cib u , mm ,Y e m (e t )Y™'*(e k ) = + £)P*,(ej • e k ). (87) 



The coefficients 6«< mm < are easy to read off and one obtains 

i>Ui mm' — ~pr 5 Hi 5mm' , 

from which one deduces 

r 1 (e 4 -e,)=^^^F,(e,-e 3 ). (89) 
This is the expression used in the main text. 
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